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ABSTRACT 

We present an evolutionary sequence of models of the photoionized disk-wind outflow around forming 
massive stars based on the Core Accretion model. The outflow is expected to be the first structure 
to be ionized by the protostar and can confine the expansion of the H II region, especially in lateral 
directions in the plane of the accretion disk. The ionizing luminosity increases as Kelvin-Helmholz 
contraction proceeds, and the H II region is formed when the stellar mass reaches ^ 10 - 20Mq 
depending on the initial cloud core properties. Although some part of outer disk surface remains 
neutral due to shielding by the inner disk and the disk wind, almost the whole of the outflow is 
ionized in lO^-lO^yr after initial H II region formation. Having calculated the extent and temperature 
structure of the H II region within the immediate protostellar environment, we then make predictions 
for the strength of its free-free continuum and recombination line emission. The free-free radio emission 
from the ionized outflow has a flux density of ~ (20-200) x (i//10GHz)^’mJy for a source at a distance 
of 1 kpc with a spectral index p ~ 0.4- 0.7, and the apparent size is typically ^500 AU at 10 GHz. 
The H40a line profile has a width of about 100 km s“^. These properties of our model are consistent 
with observed radio winds and jets around forming massive protostars. 

Keywords: stars: formation, evolution - accretion, accretion disks. 


1. INTRODUGTION 

The formation mechanism of massive stars remains 
much debated with Core Accretion, Competitive Accre- 
tion and even P rotostellar Collisions still discussed (see 
iTan et riI1l20l4 for a recent review). 

One of the key differences between high and low-mass 
protostars is that the former are expected to complete 
the later stages of their accretion while fusing hydrogen 
on the main sequence. Their photospheres should be hot 
and strong sources of Lyman continnum photons with 
energies > 13.6 eV, which can photoionize H to create 
an H II region. These protostellar H II regions are im¬ 
portant for at least two main reasons. First, by their 
long wavelength radio free-free and Hydrogen recombina¬ 
tion line (HRL) emission they provide diagnostics of the 
massive star formation process. Second, the ionized gas 
reaches temperatures of ^ 10^ K and its thermal pressure 
could be an important feedback process that regulates 
massive star formation, either by disrupting the wider 
scale accretion reservoir (i.e., core envelope or Bondi- 
Hovle competitive accretion region: e.g., iDale et aD 
I2005t [Peters et al.ll2010ll2011ll or bv photoevaporation of 

the y c retion disk (iHollenbach et al.lll994lYorke fc Weld 

19961 iRichline & Yorkel 119971 iMcKee & Tad 120081: 
l Tanakaetal. ll 2ni.€ - - - 

Our goal in this study is to present model calculations 
of the structure of the earliest-stage protostellar H II 
regions in the contex t of the Turbulent Core Model of 
iMcKee fc TanI (|2003l l (hereafter, MT03). Early work in 


this direction was carried out by iTan fc McKed (j2003ll , 
who described these structures as “Outflow-Confined” 
H II regions, since the magneto-centrifugally-launched 
protostellar outflow (either disk-wind, X-wind or both) 
w ill be t he firs t barrier to the ionizing photons. 

iKetol (120071) studied the structure and evolution of 
H II regions that are confined by the infalling flow 
onto the protostellar disk. However, their model did 
not include the effect of the MHD-driven outflow. If 
massive stars form in a scaled-up, but similar manner 
to low-mass stars, then we expect such outflows to al¬ 
ways be present around massive protostars during their 
main accretion phase. Indeed, larger-scale collimated 
outflows are oft en observed to be d r iven from massive 
protostars fe.g.. iBeuther et al.ll200l : iGarav et al.l l2007l 
iLopez-Sepulcre et al.ll2009l) 

A series of massive protostellar evolution calculations 
to predict continuum radiative transfer, especially in¬ 
frared to sub-mm morphologies and spect ral energy dis¬ 
tribut io ns (SEDs), ha s been carried out bviZhang fc TanI 
(|201lD . iZhang et al.l (|2013D . and I Zhang et al.l ( 201^ 
(hereafter, ZTH14). Here we will calculate the ioniza¬ 
tion structure and associated radio emission predicted 
by these protostellar evolutionary sequences. 

Study of radio emission f rom massive protostars has 
a long history (see, e.g., Wood fc Churchwelll 119891 : 

1201 Jb Ultra-compact (< 


iHoare et al.l l2007t ITan et al 


0.1 pc) and, especially, hyper-compact (< 0.01 pc) H II 
regions may include sources that are in the protostel¬ 
lar phase. Elongated radio continuum sources, e.g., ra- 
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dip “jets”, have been seen in a number of sources (e.g., 
iGuzman et al.l l2012f) . The eventual goal of our work 
will be to compare the developed models with such ob¬ 
servations, however, detailed comparison with individ¬ 
ual sources, requiring tailored protostellar models, is de¬ 
ferred to future works in this series. 

This paper is organized as follows. In ^we describe 
our methods: in particular, we briefly review the proto¬ 
stellar evolution models, describe calculation of the pho- 
toionized structure and the radiative transfer of free-free 
continuum and HRL emission. Next, in ^ we present 
our results: especially, we explore the evolution of the 
photoionized regions for various initial core conditions 
and predict their observational features. In 21 dis¬ 
cuss the implications of our results, especially summa¬ 
rizing the evolution of the outflow confined H II regions, 
and making an initial comparison of our models with the 
general characteristics of observed HC/UC H II regions 
and radio winds/jets. We conclude in 21 

2. METHODS 

2.1. Evolution of protostars, disks, envelopes and 
outflows 

We calculate the photoionization structure and present 
synthetic observations in the context of the Turbulent 
Core Accretion Model of MT03. We use a series of mod¬ 
els of protostellar evolution and surrounding gas struc¬ 
ture constructed by ZTH14, which self-consistently in¬ 
clude the evolution of the protostar and disk formed from 
a rotating parent core embedded in a particular clump 
environment, and mass accretion regulated by disk-wind 
outflows. We note that these models are highly idealized 
in that they assume axisymmetry about a single proto¬ 
star, while real objects will be more irregular, especially 
if the initial core is highly turbulent, and may be form¬ 
ing in a cluster environment with other protostars in the 
vicinity. Even if there is a relatively ordered and isolated 
core that is collapsing to a central disk, this collapse may 
lead to binary or low-order multiple star formation. Still, 
by first presenting the simplest case of formation of an 
isolated massive protostar, we can then see how relevant 
such models are as approximations to real systems. 

Here we briefly review the evolution of the protostar 
and its surrounding gas structures. The initial core 
is assumed to be spherical and in near virial equilib¬ 
rium via support by turbulence and/or magnetic fields. 
The evolutionary track of a core is determined by three 
environmental initial conditions: core mass Me, mass 
surface density of ambient clump Eci (which sets sur¬ 
face pressure of the core), and ratio of the core’s ini¬ 
tial rotational to gravitational energy fc- We fix the 
rotational parameter fc = 0.02 as a typical value 
from observations o f lower-mass cores ((Goodman et al.l 
iLi et alllM^ IPalau et al.l[20T^ . The core den¬ 
sity structure is described by a power law in spheri¬ 
cal radius, p oc r~^p. Values of kp ~ 1.3-1.6 are 
inferred from observat ions of dense core s in Infrared 
Dark Clouds (IRDCs) ((Butler fc Tanl[2Ml : IButler et al.l 
1201411 . We adopt kp = 1.5 as a fiducial value, fol¬ 
lowing MT03. Then, the radius of a core is Rc = 
1.2 X 1O4(Mc/6OM0)1/2(Eci/1 g cm-2)-i/2AU. Since 
the core is in pressure equilibrium with the surrounding 
clump, it is more compact for higher Ed. 
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Figure 1. Evolutionary sequence of density, tt-h, and veloc¬ 
ity fields around protostars of mass m* = 8, 12, 16, and 24Mq 
forming from a core with initial properties Me = 60Mq, Sd = 
1 g cm~^, and = 0.02 (model AOS, A12, A16, and A24). Three 
spatial scales, successively zooming out from the protostar, are 
shown from left to right. In each panel, the protostar is at the lower 
left corner, the horizontal axis lies on the disk midplane and the 
vertical axis along the outflow axis. The velocity fields are shown 
by arrows (purple for the outflow; light-blue for the infalling enve¬ 
lope; note the different velocity scales). The initial core radius, i.e., 
its boundary with the clump, and the boundaries of disk, envelope 
and outflow are indicated by white lines. 

Table [1] presents the parameter sets of (Me, Ed) for 
which we study the ionization structures. These are the 
same as studied by ZTH14. We refer to model group A 
as a fiducial group with Me = 6 QMq and Ed = 1 g cm“^. 
We also investigate high- and low-Ed cases, model groups 
A1 and Ah, and Me = I 2 OM 0 and 24OM0 cases, as model 
groups B and C, respectively. 

The core undergoes insid e-out collapse (IsbJ [I977I: 
iMcLaughlin k. Puclritzl[TM^ wi th the e ffect o f rotation 
described with the solution by lUlricbl (1197611 (see also 
iTan &: McKeel 120041) . Given the high accretion rates 
associated with collapse of dense cores, the protostel¬ 
lar disk is expected to be massive and self-gravitating 
such that angular momentum may be transported by 
torques due to, e.g., spiral arm structures. Gonsidering 
this efheient angular momentum transfer, the mass ra¬ 
tio between disk and protostar is set to be a constant 
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Table 1 

Model lists 


Models 

Me (Mq) 

Eel (gem 

m* {Mq) 

r, {Rq) 

T*,,ec (xlO^K) 

J^*,acc (X1O'‘L0) 

‘S'*,acc ( ® 

A 12 

60 

1 

12 

12.9 

2.33 

4.4 

46.5 

AI 6 



16 

6.4 

3.67 

6.6 

48.3 

A24 



24 

6.4 

3.90 

8.5 

48.6 

SIO 8 

60 

0.316 

8 

14.0 

1.59 

1.1 

43.9 

A112 



12 

4.7 

3.27 

2.3 

47.5 

A116 



16 

5.1 

3.43 

3.2 

47.8 

Ahl2 

60 

3.16 

12 

50.2 

1.21 

4.9 

42.5 

Ahl6 



16 

26.1 

2.05 

10.9 

46.1 

Ah24 



24 

9.0 

4.05 

19.6 

48.9 

BT 2 

120 

1 

12 

17.1 

TOi 

4.6 

45:7 

B16 



16 

8.4 

3.43 

8.7 

48.3 

B24 



24 

6.6 

4.20 

12.1 

48.8 

B32 



32 

7.6 

4.36 

18.9 

49.0 

CT 2 

240 

1 

12 

21.6 

1.80 

4.4 

45.2 

C16 



16 

11.1 

3.04 

9.6 

48.1 

C24 



24 

6.8 

4.30 

14.3 

48.9 

C32 



32 

7.7 

4.50 

22.1 

49.1 

C48 



48 

9.6 

4.76 

43.1 

49.4 

C64 



64 

11.5 

4.91 

68.4 

49.6 


fd = md/m^ = 1/3 fe.g.. iKratter et all I2008D . The 
disk st ructure is described by a st andard ct-viscosity disk 
model (|Shakura fc Sunva^fl973[ l. but also including the 
effects of infall and outflow (conservations of mass and 
angular momentum). The temperatures at the midplane 
are evaluated by vertical energy balance, and the vertical 
structure of the disk is determined by this temperature. 
The disk surface is set by assuming continuity with the 
density structure of the out flow. More d e tails a bout the 
disk model are described in iZhang et al.l ()2ni3[ l. 

The mass accretion rate onto the forming protostar rh* 
is 

Me / Sol 

60Mq j yl g cm“^ 

where e*(t) and e*(t) are the instantaneous and aver¬ 
aged star formation efficiencies. The instantaneous star 
formation efficiency e* is the ratio of the accretion rate 
onto the star to the idealized spherical infall rate from 
the initial cloud core in the absence of feedback: a cer¬ 
tain percentage of infalling material is swept up by the 
disk wind (see next paragraph) and a fraction, fd, of the 
accreted mass remains in the disk during the protostellar 
phase. The averaged star formation efficiency e* is the 
ratio of the stellar mass to the idealized collapsed mass 
of the initial cloud core at the same stage of infall in the 
absence of feedback. The accretion rate is high for high 
clump mass surface densities, since high Ed results in a 
more compact core and thus shorter collapse time. Also, 
this formula indicates that in the no-feedback limit of 
e* = e* = 1 , the accretion rate increases with time, i.e., 
as (see MT03). 

The magnetically driven bipolar outflows sweep up 
part of the core and thus help set the star formation 
efficiency e*. The density and velocity distributions of 
the outflow are described by a semi-analytic disk wind 
solution that is modified from the centrifu gally driven 
out flow model of jBlandford fc Pavii^ (|1982ll (for details, 
see iZhang et al.l I2013L and ZTH14b The mass loading 
rate of the wind is proportional to the mass accretion 


Meyr \ (1) 



rate, with a ratio f,,, = 0-1 assumed as a typical value 
(|Konigl fc Pudritj 120001) . The velocity is as high as 

1000 km s“^ near the rotation axis and then smoothly 
decreases at lower latitudes. This outflow velocity struc¬ 
ture represents features of both a collimated fast jet and 
a wider-angle, slower wind. The opening angle of the 
o utflow cavity ^w.bhc i s estim ated following the method 
of iMatzner &: McKed (|2000D : if the outflow momentum 
is strong enough to accelerate the core material to its 
escape speed, the outflow extends in that direction. The 
outflow momentum, p^, is evaluated by integrating the 
momentum rate over time. The momentum injection 
rate is typically Pw — (0.15 - 3) x ui^^vk* in our model, 
where vk* — 44 O(m,/lOM 0 )^/^(r*/lOi? 0 )“^/^ kms“^ is 
the Keplerian speed at the stellar radius. This process 
sets the star formation efficiency e* (and also e*). For 
more details of the outflow mod el and the evaluati on of 
the star formation efficiency, see IZhang et al.l (I2013D and 
ZTH14. 

Figure [T] shows an example of the evolution of the den¬ 
sity and velocity structures in the fiducial case. The 
gradual opening-up of the outflow can be seen at the 
scale of 2 X 10"^ AU with typical outflow velocity of 100- 
500 kms“^ which is high near the axis. This leads to 
a decrease in the instantaneous star formation efficiency. 
At the end of accretion in this model, the stellar mass 
reaches 26 Mq and the final averaged star formation ef¬ 
ficiency is e*/ = 26 M 0 / 6 OM 0 = 0.43. For the different 
models, the averaged final star formation efficiencies are 
similar, with a range of e*/ = 0.34 - 0.52. While we ex¬ 
pect the momentum flux from protostellar outflows to 
be the dominant source of feedback, we note that addi¬ 
tional feedback processes due to ionization, (line-driven) 
stellar winds and radiation pressure are not yet included 
in our modeling, so the above efficiencies may be some¬ 
what overestimated. The growth of the disk radius can 
be seen in the panels showing the scale of 100 AU. Note, 
partly for the reason of being near the end of its accre¬ 
tion phase and having a relatively low accretion rate, the 
disk in the 24 Mq case appears relatively thin compared 
to earlier stages (ZTH14). On the 100 AU scales, we also 
see the narrow, inner, on-axis cavity of the disk wind, 
which is set to have zero (i.e., negligible) density. This 
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Figure 2. Evolutionary sequences of r«, L*.acc, Tt^acc, and S».acc 
versus protostellar mass (top to bottom) for protostars forming 
from Me = 6 OM 0 and /3c = 0.02. The black solid line shows the 
case with Ed = 1 g cm“^, while blue dashed and red long-dashed 
lines show the Ed = 0.316 and 3.16 g cm“^ cases, respectively. 
The ZAMS properties are shown by the dotted lines. 

cavity is too thin to see in the larger scale panels, but it 
is resolved by the logarithmic grid spacing. 

The evolution of the protostar depends on its ac¬ 
cretion history. ZTH14 calculated the protostellar 
properties with a detailed stel lar evolution model b y 
iHosokawa fc Omukail (|2009D and iHosokawa et ^ (|2010D . 
self-consistently adapted to the Turbulent Core Model 
accretion rates. When the flow reaches the stellar sur¬ 
face, the accretion energy, Lacc = Gm*m*/(2r*), is re¬ 
leased there. Following ZTH14, we treat this accretion 
luminosity and the internal stellar luminosity as being 
emitted isotropically with a single effective temperature, 
Tlf.acc, i-e., T*,acc = i* + Face = 47rr^crT^f acc• We will re¬ 
fer to r*_acc as the stellar effective temperature. ZTH14 
also included the effects of the accretion energy released 


at the disk surfaces in the dust temperature calculation. 
However, in this study, we neglect the ionizing photons 
from the disk surface, since its surface temperature is 
not high enough to emit significant amounts of ionizing 
flux. If the protostar rotates at almost its breakup ve¬ 
locity, it distorts the structure of the star. This effect 
leads to nonuniform surface temperatures: the polar re¬ 
gion is hotter than the equatorial region (|Mevnet et al.l 
I 2 OIOD . However, the degree of rotation of massive proto¬ 
stars is quite uncertain, and thus for simplicity we adopt 
a spherical model with a single surface temperature as a 
first approximation. 

Figure [5] shows the evolution of protostellar proper¬ 
ties with protostellar mass: the stellar radius, luminos¬ 
ity and effective temperature were calculated by ZTH14, 
while the evaluation of ionizing flux will be described 
in 112.2.11 In the fiducial case with Me = 60 Mq and 
Sci = 1 g cm“^, the stellar radius reaches a maximum 
at TO* ~ 5 Mq when the Kelvin-Helmholz time be¬ 
comes as short as the accretion time, after which the 
protostar is able to begin c ontraction towards the zero - 
age main sequence (ZAMS) (IHosokawa fc Omukai2009l) . 
Due to this variation in radius, the effective temperature 
increases strongly for to* > 5Mq. Comparing to the 
fiducial case, the mass at which the stellar radius has a 
maximum is larger in the higher Ed case and smaller in 
the lower Sd case. This is because a higher Ed causes 
higher accretion rates (eq. [T]), which make the accretion 
timescale shorter and also lead to increased swelling of 
the protostar. 

2 .2. Photoionization 
2.2.1. Ionizing photon luminosity 

The total stellar H-ionizing photon luminosity S'*,acc is 
evaluated as 

= ( 2 ) 

where v is the frequency, vq is the Lyman limit frequency, 
and Li/,»,acc is the luminosity per frequency interval (i.e., 
the stellar spectrum). T he stellar atmosphere model of 
iCastelli fc Kuruci (I2004D is adopted to estimate Ljy,*,acc, 
which covers wide ranges of effective temperature and 
surface gravity. The variation of ionizing luminosity is 
sensitive to the effective temperature and varies by or¬ 
ders of magnitude (bottom panel of Figure [5]). In the 
fiducial case of group A, the ionizing luminosity is as 
small as 3 X 10^^ s“^ at 4.7 Mq, and quickly increases 
to 10^® s“^ or larger in the Kelvin-Helmholz contraction 
phase. Even at the same mass of 12 Mq, the ionizing 
luminosity of model Ahl2 is five orders of magnitude 
smaller than that of model A112. Therefore, the accre¬ 
tion history and the detailed stellar evolution calculation 
are essential to study the ionization structure. We note 
that the obtained ionizing photon lumi nosities are some¬ 
what lower than those in Figure 2 of iTan et al.l (120141 1 
due to improved calculation of the effects of line blanket¬ 
ing. 

2.2.2. Transfer of ionizing photons 

We calculate the extent of the ionized regions us- 
ing our previously developed radiative transfer code 
(| Tanaka et al.l 1201^ . which allows a treatment of both 
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direct and diffuse ionizing radiation fields. We consider 
only the ionization of Hydrogen, ignoring Helium and 
other elements. If the stellar effective temperature is 
high, T*_acc > 40, OOOK, we would underestimate the elec¬ 
tron density by about 10% due to singly ionized Helium, 
which would then enhance free-free emission by about 
20 %. 

The gray approximation for the photoionization calcu¬ 
lation is adopted, where the frequency of ionizing pho¬ 
tons is represented by the mean value. The equation of 
radiative transfer along a ray is 

dl T 2 2 

— = -xinacrul -I- aiXiiJiYi— — 
dS TTT 

-n-HO-d,a4^ + n-HO-d,s(^-4), (3) 


where nn is the number density of hydrogen nuclei, / is 
the irradiance intensity, J = J Idfl/in is the mean in¬ 
tensity, {hiy) is the mean energy of ionizing photons, cth 
is the mean photoionization cross-section of the Hydro¬ 
gen atom, ai is the radiative recombination coefficient 
to the ground state, xi and a;ii are neutral (atomic) and 
ionized fractions of H, crd,a and (Td,s are absorption and 
scattering cross sections of dust grains per H nucleon. 
The first term of the right-hand side represents the con¬ 
sumption of ionizing photons by neutral Hydrogen, the 
second term describes the re-emission of ionizing pho¬ 
tons by recombinations directly to the ground state, and 
the third and fourth terms correspond to absorption and 
scattering by dust grains, respectively. The local ion¬ 
ization fraction of H, xu, is calculated from the balance 
between photoionization, recombination and advection. 


V • (nHa;iiv) 


inxinnauJ 

{hv) 


— aAXiin^, 


(4) 


where v is the velocity of the disk wind gas, aA is the 
recombination coefficient for all levels (so-called case A). 
The advection term on the left-hand side of this equation 
can be especially important near the ionization bound¬ 
ary. The recombination rates, ai and a a, are functions 
of temperature of the ionized gas. The mean values of the 
photon energy and Hydrogen cross-section are evaluated 
from the stellar spectra for each model, 


(hi/) 


1 

'*,acc 


1 


CTH = 


S 


*,acc 



dv, 


(5) 

( 6 ) 


where cth,^ is the frequency dependent cross-section of 
the Hy d rogeii atom. We adopt the fitting formulas from 
iDraind (120111 1 for the recombination rates ai and aA, 
and the Hydrogen cross-section cth,!^. The absorption 
and scattering properties of dust in these H II regions 
are quite uncertain. We adopt CTd,a and CTd,s equal to 
zero in regions where the temperature history has ex¬ 
ceeded the dust destruction temperature (1,400 K) in 
the thermal radiation calculation by ZTH14 (i.e., there 
is a dust-free inner region of the disk wind), and else¬ 
where they are fixed at 10“^^ c m^ from a typical value 
of the diffuse interstellar medium (|Weingartner fc Drainel 
1200111 . The effect of variation of dust grain properties, 
i.e., sizes, abundances, compositions, is deferred to a fu¬ 
ture paper. 


We conduct axisymmetric two dimensional calcula¬ 
tions to solve equations and 0, using our devel¬ 
oped ra y-tracing radiative t ransfer code. In the previous 
work of l Tanaka et al.l (1201311 , the code solved the transfer 
equation by the method of short characteristics, which 
involves only the cells nex t to the point under consid¬ 
eration (|Stone et al.l 1199211 . Here, we update the code 
implementing the method of long characteristics. With 
this method, the transfer problem is solved along the 
rays from the domain boundary to each point, which 
leads to less numerical diffusion and treats the diffuse 
radiation more accurately than the short characteristics 
method. The computational domain is a cylindrical re¬ 
gion with radial and vertical coordinates R < 20 ,000 AU 
and Z < 20,000 AU, assuming equatorial plane symme¬ 
try. The ionizing photon source is a sphere located at 
{R, Z) = (0,0) with a radius of r,, which represents the 
accreting star. We denote the distance from the center 
of the star as r = + Z"^. The spatial grids for radial 

and vertical directions are set to be spaced logarithmi¬ 
cally since the inner region has finer structures. The en¬ 
tire computational domain is resolved with a 256 x 256 
grid. The radiation intensity / is a function of location 
and also direction: I = I{R,Z,9,(j)) where 9 and (p are 
the polar and azimuthal angles, respectively. The radia¬ 
tion emitted from the vicinity of the star has very fine an¬ 
gular structure when it reaches the outer boundary, i.e., 
69 ^ S(j) ^ r*/ri„ax 10“®. To resolve these fine struc¬ 
tures while still saving computational costs, the angular 
resolutions of 9 and (p are configured to vary depending 
on location {R, Z). In the 9 direction, the radiation 
arriving from directions with 9 ~ tan“^(i?/Z), which 
comes from the vicinity of the star, is resolved with a 
high resolution of A9 <C r^/r, while the radiation arriv¬ 
ing from 9 ~ tan“^(i?/Z)-|-7r/2 is resolved by A9 ~ 10“^. 
In the (p direction, based on the tangential plane method 
(|Stone et al.lll992l ). the radiation emitted from near the 
axis and from the opposite direction, ~ 0 and tt, is 
resolved with high resolution of Acp <C r*/r, while the 
radiation from (p ~ tt/2 is resolved by A(p ~ 10“^. We 
also conduct resolution tests for the fiducial case with 
grids of 128 x 128 and 330 x 330. The obtained ionized 
structures are consistent at these resolutions and free- 
free fluxes change by less than 6% at radio frequencies 
from 0.3-300GHz. 

In this hrst work we study the photoionization of the 
disk wind, and do not allow photoionization of gas in the 
disk or infall envelope components of the model. Such 
photoionization would enhance the mass loss rate that 
is being injected into the disk wind, i.e., by photoevapo¬ 
ration, and would thus change its density structure and 
also the star formation efficiency. This feedback of pho¬ 
toevaporation will be discussed in a future paper. 

2.2.3. Temperature of the ionized gas 

The temperature of the ionized gas is important to 
evaluate recombination rates and to study free-free emis¬ 
sion for observational diagnostics. The photoionized gas 
temperature is regulated to around 10, 000 K by the bal¬ 
ance of heating and cooling processes. The dominant 
heating process is photoionization and the main cool¬ 
ing processes are radiative cooling of collisionally excited 
metal lines, recombination of H, and free-free emission. 
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Figure 3. Top: Radial structure of an H II region with an uni¬ 
form density of njj = 10® cm“® formed by the radiation of the 
stellar model of A24. The radial structures of temperature, T, 
(solid) and ionization degree of H, xn, (dashed) are shown. The 
volume averaged temperature of Tave = 10, 225 K is also indicated 
(dotted). Bottom: Ionized gas temperature as a function of density, 
T(nH), which is set equal to the volume-averaged photoionized gas 
temperature, Tave(>^H)- Results are shown for stellar models A12, 
A16, and A24. 


It is computationally expensive to take into account all 
heating and cooling processes at every step of the radia¬ 
tive transfer calculation, and thus here we introduce an 
approximate treatment. The heating rate by photoion¬ 
ization depends on the stellar spectrum. Also, heating 
and cooling rates are both mostly sensitive to the den¬ 
sity. Therefore, we treat the photoionized gas temper¬ 
ature, for each stellar model, as a function of the local 
density, i.e., T{npi). We evaluate T(nH) from the volume 
average temperature of static H II regions with uniform 
densities in the ran ge nu = 10 - 10^^ cin ~^ using version 
13.03 of CLOUDY l)Ferland et al.l 1201^ . This approx¬ 
imate approach is appropriate since a uniform density 
H II region for each stellar model is almost isothermal 
as shown below. It should be noted that the ionized gas 
temperature is also sensitive to the metallicity of the gas, 
and we assume solar metallicity in this study. 

The top panel of Figure |3] shows the radial structure of 
an H II region with uniform density of nn = 10® cm“® 
formed by the stellar radiation of the stellar model A24. 
It can be seen that the H II region has a sharp edge at 
960AU where the ionization degree drops from near unity 
to near zero. The temperature is almost isothermal in the 
ionized region. We evaluate the typical photoionized gas 
temperature at nn = 10® cm“® as being 10, 225K from 
the volume average of this temperature for this stellar 
model. 

In this way, we obtain photoionized gas temperatures 
over the density range of tih = 10 - 10^^ cm“® for each 
stellar model. The bottom panel of Figure |3] shows the 


photoionized gas temperature as functions of density, 
T(nH) = T'ave(nH), for the stellar models A12, A16, and 
A24. The temperatures of all models are around 5, 000- 
15,000K for orders of magnitude of density variation. 
The temperature is typically higher at higher density. 
This is because the cooling is suppressed if the density 
exceeds the critical density of some of the cooling lines. 
We adopt these averaged temperatures to set T{nH) as 
the local temperature of the ionized gas in our radiative 
transfer calculations. 

Since the adopted temperature is evaluated based on 
static H II region calculations, adiabatic cooling has not 
been included in this study. The adiabatic cooling rate, 
Aad — 2.5V • (npikBTv), is proportional to ~ rtBlr. On 
the other hand, the radiative cooling rates in ionized re¬ 
gions are almost proportional to the square of density. 
Since the density gradient in the disk wind is as steep 
as dlogrin/dlogr ~ —1.5 to —2, adiabatic cooling may 
not be negligible in low density regions far from the star. 
We compare the adiabatic cooling rate to the cooling 
rate evaluated by the static calculations. In the fiducial 
model of A16, adiabatic cooling would be > 20% of total 
cooling rate when nn < 10^ cm“® and r > 1000 AU. As 
we will see in H3.2.11 the free-free flux at frequencies lower 
than 1 GHz is mainly provided from r > lOOOAU. There¬ 
fore, we note that the SEDs obtained by our simplified 
modeling would be overestimated somewhat at < 1 GHz. 


2.3. Observables 
2.3.1. Free-free emission 

Using the obtained structures of the photoionized out¬ 
flow, we estimate the free-free radio continuum emission. 
The equation of transfer for free-free emission is 


{I.,s - B ,), (7) 

as 

where and are the intensity and opacity of free- 
free emission, and B^{T) is the Planck function. The 
opacity due to free-free absorption is given by 


4 /27r 




1/2 


^H^HII 




m, 


1 — exp {—hv/ksT) 




( 8 ) 


where e and rrie are the charge and mass of the electron, 
h is the Planck constant, c is the speed of light, kg is the 
Boltzmann constant, and gsi^, T) is the Gaunt factor for 
free-free emission. The observed flux density is obtained 
by integrating over source solid angle of Hs, 

7% ff = / COs(^ ^view)dll = / (9) 

dna dOs 


where 0view is the viewing inclination measured from 
the outflow axis. In the last formula, we assumed that 
the distance is much farther than the object size, i.e., 
Q — ^view ^ 20, OOOAU/1 kpc <C 1. In this study, for 
quantitative results, we assume protostars are located at 
a distance of 1 kpc, and total fluxes are obtained by 
integrating over a radius out to 10'^ AU from the cen¬ 
tral star (for all our models this region dominates the 
total flux from free-free emission). The free-free emis¬ 
sion is typically observed at ~ 0.1-100 GHz, while 
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the thermal dust emission would usually be dominant 
at higher frequencies. However, for completeness, we 
calculate free-free emission over the frequency range of 
V = 10-1-10^ GHz. 


2.3.2. Hydrogen recombination lines 

We also evaluate the HRL emission from the photoion- 
ized outflow. Considering a Hna line, i.e., the line from 
the transition n + 1 ^ n, the opacity is 


Nj^.HRL — 


^ fn + 1 
Stt 


( 10 ) 


where (j)i, is the normalized line profile = 1), 

n„ is the population number in energy level n, 
is the Einstein coefficient of spontaneous emission, and 
hvn and A„ = c/u„ are the energy and wavelength of 
the Hna transition, respectively. Here, for simplicity, 
we assume the case of local thermodynamic equilibrium 
(LTE). However we note that, in the particular case 
of the massive star MWC 349A, maser amplification of 
HRLs are observed and thus non-LTE effects can be im¬ 
portant (fo r modeling of MWC 349 A including non-LTE 
effects, see iBaez-Rubio et al.ll2013l l. 

The transfer for the total intensity, i-e., the 

sum of the free-free and HRL emission, is 


d/^,ff,HRL 

ds 


(Kj/.S + Kj^^hRl) (G,ff,HRL — Si/) , (11) 


and the total observed flux is 


Ai/,ff3RL = / 7i/,ff,HRL cos{9 — 9vie^)dV, 

J Hs 

= [ Ir,S,RRLd^. ( 12 ) 

J Oa 

Here, we assumed that the source distance is greater than 
its size as in equation (l9|). We evaluate the HRL flux as 
the excess from the continuum flux at Vn, i-e., EI^^rl = 
Fi^,s,hrl — Frn,s- Of course, observationally it is not 
possible to observe the free-free flux exactly at since 
the line exists there and we cannot divide the total flux 
into two components. However, the value of is still 
able to evaluated from the far tail of the line profile or 
from interpolating from continuum observations at other 
frequencies. Since the typical line width is Au/u„ 
vfc < 10“^ and the typical free-free emission spectral 
index in this study is 0 < p < 1, the variation of the 
continuum flux in this width is less than 0.1%. Thus, 
it is safe to assume the continuum is constant over the 
line width. In this paper we will illustrate these results 
with the H40a line at 99 GHz. The flux ratio is typically 
hrl/F^,s ^ 0.1 for the H40a line in our models. 

We assume that locally the gas has a Maxwellian ve¬ 
locity distribution, ignoring the intrinsic width of the 
line relative to the thermal and bulk-motion broadening. 
The velocity dispersion of Hydrogen is ay = \/k^T/mu, 
and the corresponding dispersion in frequency is a^, = 
{ayjc)un. The local normalized line profile pi, combin¬ 


ing the bulk gas velocity along the line of sight, v, is 


pL- 


1 

\/27r cr^ 


exp 


{(1 -I- v/c)iy - 

2al 


(13) 


The frequency v observationally represents a radial ve¬ 
locity of (1 — v/vn)c. 


2.4. Escape fraction of ionizing photons 

As we will see in next section, the ionized region, af¬ 
ter initial early expansion, becomes unconfined in the 
radial direction along the outflow axis. Thus, some frac¬ 
tion of the ionizing photons escapes from the protostel- 
lar core and can ionize the surrounding, ambient clump 
gas, which may be forming a stellar cluster. Although 
the main focus of this work is the ionized structure of 
disk winds and their observational properties, ambient 
gas ionization is also important for its effect on star clus¬ 
ter formation and production of additional radio emis¬ 
sion. Therefore, we provide estimates of the ionization 
of escaping photons and radio emission from the ionized 
clump gas. 

The escape fraction of ionizing photons depends on the 
polar angle 9, i.e., in response to the ionized disk wind 
structure. We calculate the escape fraction as 


47rr^Fr,out(^) 


*S^=(',acc(^u) 


(14) 


where E).^out is the radial energy flux of ionizing pho¬ 
tons at the outer boundary, i.e., R = 20,000 AU or 
Z = 20,000AU. We note that the escape fraction is 
insensitive to the choice of the radius at which escape is 
defined as long as it is significantly larger than the disk 
radius. Integrating /esc over the boundary, we obtain the 
total escape fraction. 


/ 


esc 


1 

47r 


fesc{0)dn. 


(15) 


3. RESULTS 

3.1. Structure and evolution of the photoionized outflow 

Eirst we report results of the fiducial case, model group 
A. Figure S] shows the evolution of the photoionized 
structure. The high temperature region with 10^ K 
is the photoionized gas, i.e., the H II region, and else¬ 
where is neutral gas (with temperatures calculated from 
the ZTH14 radiative transfer models; note, we expect 
there would be some localized changes to the ZTH14 
results for the neutral gas due to the presence of the 
H II region and its associated photodissociation region 
(PDR), but we ignore those here—the PDR region will 
be considered in a future paper). Note that the white 
regions seen along the rotation axes in the 100 AU-scale 
panels are disk wind inner cavities, where the density is 
negligible (formally zero in these models since we do not 
include a line-driven wind from the stellar surface) and 
thus where ionizing photons can travel freely. The ion¬ 
ized gas temperature is slightly higher in the inner region 
of the disk wind due to the higher density (see Fig. [S]). 
When m* < SMq, the ionizing photon luminosity is too 
low to ionize the disk wind. As the ionizing luminosity 
increases with mass, the photoionized region is formed. 
When TO* = 12 Mq, the photoionized region is confined 
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Figure 4. Evolutionary sequences of the photoionized outflows 
of model group A. The color scale shows the temperature of the 
photoionized and neutral regions. The white regions seen close to 
the vertical axes in the 100 AU-scale panels (left column) are disk 
wind inner cavities, where the gas density is negligible (formally 
zero) and the ionizing photon can travel freely. As in Figure^ the 
core scale and the boundaries between disk, core envelope, outflow 
and clump are shown by white lines. 

to intermediate zenith angles of 10° ^ 6* < 30°. This 
is because densities are higher near the outflow axis at 
0 < 10° and near its launching zone for 0 > 30°. Since 
the radial density gradient in the outflow is typically 
steeper than the photoionized region is not con¬ 

fined in the radial direction and breaks out. The opening 
angle of the photoionized region increases as the proto- 
stellar mass grows. The ionized region eventually reaches 
the flared disk surfaces, however some parts of the disk 
wind remain neutral due to the shielding by the inner 
disk and disk wind even though the diffuse radiation is 
considered. Therefore, photoevaporation from the disk 
surface is suppressed, while photoevaporation from the 
envelope may still be efficient. 

The ionizing radiation has two components: one is 
the direct stellar radiation and the other is the diffuse 
radiation from recombinations to the ground state and 


le2AU leSAU 2e4AU 



Jdiff/(J.+Jdifif) 

Figure 5. The radiation fields for model group A. The color 
scale shows the fraction of diffuse radiation in the ionized region 
(the white region is neutral). Black lines indicate the disk/envelope 
boundaries with the disk wind. 


log R [AU] log R [AU] 

Figure 6. The structures of density uh (left) and ionization 
degree xu (right) for model A16 in logarithmic scale. The white 
region in the right panel is the inner hole where the ionization 
degree cannot be defined due to the zero density. However, the 
ionizing photons travels this inner hole freely. Note that the stellar 
radius is r* = 6ARq = 0.03 AU for this model. 

dust scattering. Figure O shows the fraction of diffuse 
radiation, Jdis/{J* + ./diff), at each evolutionary stage. 
The direct radiation dominates over most of the outflow. 
However, in some locations, due to the direct radiation 
being shielded by the inner disk at -C 10^ AU, the diffuse 
radiation becomes more important, i.e., near the ionized- 
neutral boundar y, e.g., as seen in model A 24. In this way, 
as suggested bv iHollenbach et all (|1994ll . the diffuse ra¬ 
diation is important to determine the structure of the 
H II region and rates of photoevaporation from the neu- 
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tral components. For these reasons, accurate treatment 
of the radiative transfer including the diffuse component 
is necessary. 

To illustrate the structure near the rotation axis, Fig¬ 
ure [6] shows the density and ionization fraction of model 
A16 with a logarithmic scale. The disk wind inner cav¬ 
ity, with zero density inside its innermost streamline, is 
seen as a very collimated structure with 9 < 2 x 10“^ at 
Z = 20,000 AU (the black region in left panel and the 
white region in right panel). The density is high along 
the inner edge of the outflow. Due to the shielding by 
this high density inner wall, the gas in the region at high 
latitude with R ~ 10-100 AU and Z > 1000 AU remains 
neutral. This shielding is the same mechanism seen in 
model A12 in which the neutral region is at 0 < 10°. The 
strip of ionized gas entering this neutral region is caused 
by high velocity advection near the axis (eq. H]). These 
thin structures near the axis are sensitive to our adop¬ 
tion of an idealized axisymmetric model. In reality, the 
outflow winds are not expected to be perfectly axisym¬ 
metric. For example, t he jet in IRAS 16547-42 47 appears 
to undergo precession (jRodrfguez et al.ir2008ll . However, 
the idealized structures in our model illustrate the wind 
shielding effect and also demonstrate the accuracy of our 
numerical calculation even at such fine angular scales. 

Next we show results with different initial cores with 
lower and higher Ed. The evolution of density and pho¬ 
toionization (as revealed by gas temperature, T) struc¬ 
tures in groups A1 and Ah are shown in Figures [7] and 
[HI In group Al, the photoionized region breaks out when 
8Mq < m* < I2M0, as in the fiducial case. On the other 
hand in group Ah, the photoionized region is tightly con¬ 
fined to a thin layer at the inner wind wall even at I6M0 
and only expands later. This is mostly because higher 
Ed leads to higher accretion rates (eq. |T]), which result 
in a higher mass at which the ionizing luminosity shows 
its dramatic increase due to Kelvin-Helmholz contraction 
(Fig. [5]). We discuss the breakout of outflow-confined HII 
regions in 114.11 

Figures 111 an do] show the evolution of density and pho¬ 
toionization (as revealed by gas temperature, T) struc¬ 
tures in higher Me models of groups B and C. In both 
cases, the photoionized region breaks out when m* > 
12 M0. This is similar to the high Ed case, as the higher- 
mass core leads to delayed (in m*) contraction via higher 
accretion rates. If we ignore feedback effects, the accre¬ 
tion rate is proportional to at the same stellar mass 

(eq. [IJ ■ In this way, the mass at which the photoionized 
region emerges varies from ~ 10-15 Mq depending on 
its parent cloud core. 

3.2. Free-free continuum emission 

We conduct the radiative transfer calculation of free- 
free emission from the obtained photoionized structures. 
The total flux, spectrum, and size of the free-free emis¬ 
sion of our model are similar to those of the radio winds 
and jets associating with massive star forming regions as 
we will see in 94.21 However, our ionized region model has 
a limited cylindrical volume (' 92.2.2|1 and the contribu¬ 
tion from external ionized regions is not included. Thus, 
it should be kept in mind that the free-free flux would be 
enhanced by the emission from the ionized clump if it is 
on the line of sight, i.e., 0view < ^w.esc- 


3.2.1. Spectral Energy Distributions 

For all SED calculations, a distance of I kpc and an 
aperture diameter of 20" (20,000 AU) are assumed. Fig¬ 
ure HD shows full SEDs including thermal dust and free- 
free emission for model A16 at viewing inclinations of 60° 
and 30° of the line of sight to the rotation/outflow axis. 
The flux from thermal dust emission was calculated by 
ZTH14. The free-free emission dominates over the ther¬ 
mal dust emission at radio frequencies of < 100 GHz, 
while the dust flux dominates at higher frequencies. We 
note that, at wavelengths shorter than 100 /rm, we do 
not count the free-free emission in the total flux since we 
do not treat the dust extinction in the free-free radiative 
transfer calculation (eq. [7]). 

The infrared (IR) dust flux (at relatively short wave¬ 
lengths) depends strongly on viewing angle. This is be¬ 
cause the stellar radiation is absorbed by the dense disk 
and core envelope along directions close to the equato¬ 
rial plane, and the re-emitted IR radiation, still suffering 
from high optical depths, escapes along the polar, low- 
density outflow cavity. Therefore, the observed flux is 
higher if the line of sight is passing through the outflow 
cavity (ZTHI4). This process is known as the “flash¬ 
light effect,” which is also importan t for helping to over¬ 
come radiation pressure feedback (iNakano et al. W95; 
Yorke fc BodenheimedUO^ iKrumholz et al.ll2005 , 2009 : 
Kuiner et abll201 Ol 1201 IF 

On the other hand, the radio flux is almost indepen¬ 
dent of viewing angle since it typically experiences much 
smaller optical depths than the IR radiation. We can ap¬ 
proximate the spectrum of free-free emission as being a 
power law in frequency, oc with spectral index p. 
From 3 - 30 GHz, the free-free flux from model AI6 is fit¬ 
ted as = 81(u/GHz)° "‘® mJy. The spectral index of 
the free-free emission decreases from p = 0.9 at 0.1 GHz 
to 0.2 at 100 GHz. An index smaller than two indicates 
that the photoionized disk wind is only partially optically 
thick with respect to its free-free emission. The overall 
spectral index (including thermal dust emission) has a 
minimum of 0.45 at 35 GHz, and quickly increases at 
higher frequencies approaching the dust spectrum value 
of p = 3.3. Figure [TT] also indicates that the free-free 
flux at ^ 0.1 GHz is mostly emitted from ~ 10^-10^ AU 
scales, while that at ^ 10 GHz is mostly probing ^ 10^- 
10^ AU scales. Thus, the higher frequencies trace the 
structure of the inner wind. 

Figure [12] shows SEDs for all models. The features 
of the free-free emission are similar to those discussed 
for model A16: the free-free fluxes at radio frequencies 
are almost independent of viewing angle, and are ^ 20 
- 200 mJy with the spectral indices of p ~ 0.4 - 0.7 at 
10 GHz. These common features are caused by the fact 
that the outflows of all models have similar structures, 
i.e., the amount of gas and the density gradient. If the 
photoionized region is small, the free-free emission is not 
observable, as in, e.g., model A108. Once the photoion¬ 
ized region breaks out, the free-free emission becomes 
bright, as in, e.g., model A112. Since the flux of free-free 
emission depends on the density structure of the out¬ 
flow, the mass of the outflow that is emitting free-free 
flux (within a given distance from the protostar) may be 
estimated from its free-free flux. 

With our focus in this paper on radio wavelengths. 
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Figure 7. The evolution of density (left) and photoionization (as revealed by gas temperature, T) (right) structures for model groups A1 
(Me = 60 Mq, Eel = 0.316g cm“^). 


we note again that we have not accounted for the free- 
free emission for the total flux at wavelengths less than 
100 /im. However, in some cases the free-free flux might 
make a significant contribution in this range: e.g., at 
wavelengths shorter than 10/rm of model C16 viewed at 
60°. The potential importance of free-free emission at 
shorter wavelengths will be explored in a future paper. 

3.2.2. Images 

Next, we present images of the free-free emission. We 
show both resolved images and those after convolving 
with finite beam sizes of various telescopes. Detailed 
comparisons of models with individual sources are de¬ 
ferred to future papers in this series. 

The top three panels of Figure [13] show the resolved 
free-free emission images of model A16 at the incli¬ 
nation of 60° at three different frequencies, 1.5, 10, 
and 870 GHz. These frequencies correspond to the 
L and X bands of the Karl G. Jansky Very Large 
Array (J-VLA) and band 10 of the Atacama Large 
Millimeter/submillimeter Array (ALMA). Box sizes are 
20,000 AU X 20,000 AU, 4,000 AU x 4,000 AU, and 
100 AU X 100 AU, respectively. The apparent size of 
the free-free images is larger at lower frequencies since 
free-free absorption is stronger (eq. |8|). The effec¬ 
tive radius at which half of the total flux is emitted is 
10-3(zz/10 GHz)« pc ~ 2 X 102(zz/10GHz)« AU, with 
q ~ —0.7 for 1-100 GHz. Thus these would typically 
be categorized as HG H II regions (see also M.2I) . At 
1.5 GHz and 10 GHz, the apparent shape is similar to 
an hour-glass shape with a jet, because the disk wind 


has a self-similar structure at the scale of > 1000 AU. 
The surface brightness is high in the central wind region 
and at the jet axis where the density is high, and the 
maximum brightness temperature of Tb ~ 10"^ K indi¬ 
cates optically thick conditions. At the higher frequency 
of 870 GHz, the shape is now more rounded in the inner 
region. This is because of the accretion disk with ra¬ 
dius of ^ lOOAU. The collimated on-axis jet is not seen, 
since it becomes optically thin at this frequency. For the 
total free-free flux in Figure 1111 the innermost wind re¬ 
gion dominates the on-axis jet since its area is bigger. 
However, the bright linear feature is still very interest- 
ing, since it is po t entially related to obs e rved radio jets 
(iGibb et al.l 120031 : iRodriguez et al.l [20081: iGuzman et al.l 
I2012LI2014D . 

The bottom three panels of Figure [T3| show the free- 
free images convolved with the finite beam sizes of the 
telescopes at these frequencies. The full-width at half¬ 
maximum (FWHM) beam sizes are 1.3, 0.2, and 0.01 " 
for 1.5, 10, and 870 GHz, respectively, i.e., the highest 
resolutions in the L and X bands of the J-VLA and in 
band 10 of ALMA. A distance of 1 kpc is assumed, and 
the point-spread function is approximated as a Gaussian 
function. The jet feature with high brightness tempera¬ 
ture along the axis is diluted at 1.5 GHz and 8.4 GHz, 
since it is narrower than the beam size. However, the 
linear feature is still observable at the level of T;, < 30 K 
or lOO/rJy beam“^ at both frequencies. The sensitivity 
with root-mean-square value of ^ lO^Jy beam”^ can be 
attained with an observing time of 1 hour and a band- 
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Figure 8. The evolution of density (left) and photoionization (as revealed by gas temperature, T) (right) structures for model groups 
Ah (Me = 60 Mq, Eel = 3.16g cm“^). 


width of 500 MHz by the J-VLA. The inner brightest 
region with Tf, ~ 10^ K is smaller than the beam size, 
and thus the highest brightness temperature decreases 
to 1600 K at 870 GHz. The dust emission cannot have a 
brightness temperature above the dust sublimation tem¬ 
perature of 1400 K. Therefore, even though the total 
flux is dominated by dust emission at the frequencies 
> 100 GHz (Fig. [TT|) . the free-free emission may be ob¬ 
servable even at this high frequency. Since most massive 
protostars are at somewhat larger distances, this emis¬ 
sion may need to be observed in more massive, denser 
sources. 

Figure (TT] shows resolved and convolved images at 
10 GHz as the protostar evolves, viewed at inclination 
angles of 30° and 60° for all models. From the resolved 
images (top panels of Fig. [14]), it can be seen that the 
outflow opening angle gradually increases, especially at 
the inclination of 60°, once the photoionized region is 
formed. When the line of sight passes through the out¬ 
flow cavity (dview < fi*w,esc), Hke model A16 at 30° inclina¬ 
tion, then the image becomes elliptical rather than hour¬ 
glass shaped. The outflow axis always exhibits a high 
brightness at this frequency. As stated above, the linear 
jet-like feature is observable at the level of Tf, < 30 K 
in the convolved image of model A16. In the case of 
higher-Sci or Me, the brightness temperature is higher 
since the outflow rate is higher. Thus the jet feature is 
observable even at T;, ~ 100 K in these models, such as 
models Ah24 and G32. The small spots in convolved im¬ 
ages of models Ahl6 and B12 are the sign of unresolved 


outflow-confined H II regions. 

3.3. Hydrogen recombination lines 

We perform radiative transfer of the H40a line to ob¬ 
tain simulated observations that probe the kinematics of 
the photoionized wind. Here, a distance of 1 kpc and 
aperture diameter of 20" (20,000 AU) are assumed. 

Figure [T51 shows profiles of H40q; flux from model A16. 
The line profiles are very similar for the inclination an¬ 
gles of 60° and 30°: the FWHM are 104 and 115kms“^, 
the peak fluxes are 29 and 27 mJy, the ratios of line flux 
to the free-free continuum flux are 0.16 and 0.18, respec¬ 
tively. Although we used a Gaussian line profile of ipu at 
each position for the radiative transfer calculation (eq. 

the obtained total spectrum has broader wings. This 
is because the disk wind has higher velocity components, 
> 500 km s“^, near the axis. 

Figure [T6| shows resolved and convolved maps of 
H40a of model A16 at three different velocities, 
— 125, 0, and 125kms“^, at inclinations of 60° and 30°. 
The map is convolved with a beam size of 0.1 " which is 
the highest resolution in band 3 of ALMA. For this ion¬ 
ized bipolar outflow, the blueshifted component appears 
on the north (upper) side and the redshifted component 
in south (lower) side. However, the 0kms“^ maps are 
tilted slightly to northwest-southeast direction. This is 
due to the rotation of the outflow. Rotating outflows 
have been f ound around low - mass and intermediate-mas s 
protostars (|Klaassen et al ] ICodella et afl (Ml . 

but not yet for high-mass objects. In resolved maps at an 
inclination of 60°, the line of sight is outside the cavity 

































12 


Tanaka, Tan, & Zhang 






(N 


fN 

m 




B : I 2 OM 0 , 1 g cm-2 

le2AU le3AU 2e4AU le2AU le3AU 2e4AU 


I000km 


I 


HH 




0 0.5e4 le4 1.5e4 2e4 



T[K] 


Figure 9. The evolution of density (left) and photoionization (as revealed by gas temperature, T) (right) structures for model groups B 
[Me = 120 Mq, Eel = Ig cm“^). 


and thus the rim of the outflow is brightest. On the other 
hand, in resolved maps at an inclination of 30°, the line of 
sight passes through the cavity and the maps are ellipti¬ 
cal illustrating the outflow cone. This difference becomes 
less clear in the convolved maps at ±125kms“^. How¬ 
ever, it is still visible in the 0 km s“^ maps at the level of 
Tf, ~ lOK. The sensitivity of = 2K could be attained 
with an observing time of 8 hours with a channel width 
of 6 MHz by ALMA. The effective spectral resolution 
is two times of the channel width which corresponds to 
12 MHz or 36kms“^ (the H40a HRL is at 99 GHz). We 
should mention that again a distance of 1 kpc is assumed 
here, and most massive protostars are typically at larger 
distances. Thus, this emission may need to be observed 
in more massive and denser sources. 

3.4. Escape fraction of ionizing photons 

So far we presented results for the ionization of disk 
winds. Next, we consider the ionizing photons that es¬ 
cape along the outflow cavities created by these winds. 
Since the obtained ionized regions are unconflned in the 
radial direction, some fraction of ionizing photons escape 
and are expected to ionize the ambient clump gas. This 
ionization feedback, in addition to that of the disk wind 
protostellar outflows, may affect star formation around 


the forming massive star and will also lead to additional 
radio emission. 

The top panel of Figure [T7] shows the escape fraction 
of ionizing photons to each polar angle, /esc(0), along 
with the opening angle of the disk wind cavity, dw.esc- 
In the case of model A12, the escape fraction is 0.35 at 
its maximum value and almost zero for 9 < 10° or d > 
30°. This is because the ionized region is restricted to 
10° ^ ^ 30° (Fig. 2]). As the ionized region expands in 

models A16 and A24, the angular range in which photons 
escape also expands and the maximum value of /esc(0) 
reaches about unity. However, the escape fraction is still 
< 0.5 in directions near the boundary of the core infall 
envelope (i.e., 9 > 0w,esc — 10°: see also Fig. [5]) or near 
the outflow axis {9 < 5°: see also Fig. [01) due to strong 
attenuation in these higher density regions. 

The bottom panel of Figure [T71 shows the evolution of 
the outflow cavity opening angle, dw.esc, and the total 
escape fraction of ionizing photons, f^^ci all model 
groups. It can be seen that the total escape fraction 
increases as the opening angle grows in every case. How¬ 
ever, we note that, due to Hydrogen ionization and dust 
absorption, the total escape fraction is lower than the 
fraction of solid angle extended by the disk wind cavity. 
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Figure 10. The evolution of density (left) and photoionization (as revealed by gas temperature, T) (right) structures for model groups 
C (Me = 240 Mq, Eel = Ig cm“^). 


The maximum value of in our models is 0.38 with 
^w,esc = 72°, so more than half of ionizing photons do not 
escape from the system (including the directions blocked 
by the disk and infalling core envelope). 

The total ionizing photon rate escaping from the pro- 
tostellar core is /esc'S'*,acc. Most of these photons will 
be able to ionize the ambient clump or larger-scale cloud 
gas, especially the gas in front of the head of the out¬ 
flow bowshock. Our model does not explicitly include 
the structure of the outflow bowshock or the ambient 
clump gas. Thus, we do not study details of the geom¬ 
etry or kinematics of the ionized clump. Instead, as a 


simple hrst approximation, we use a Stromgren sphere 
analysis to make some simple estimations. In this case, 
the approximate size scale of the ionized region created 
by escaping photons is 


7?st,ci 


3/, 


esc^*,acc 


1/3 




4«b 


(16) 


where nn.ci is the density of the ambient clump, and 
Ob = CIA — cii is the case B recombination rate 
(this estimate ignores absorption of photons by dust). 
The mean clump density is given by uh.cI = 2.0 x 
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Figure 11. SEDs for model A16 {Me = 60 Mq, Ed = 1 gem and = 16 Mq) at inclination angles of 60° (left) and 30° (right). 
Two emission components, free-free (thick blue) and thermal dust (thin black), are shown in each panel. The blue dotted lines show the 
contributions by the regions at 1,10,10^, 10^ and 10^ AU from the central star. A distance of 1 kpc is assumed for the flux normalizations. 


lOySci/l g cm-2)3/2(Mci/4000 cm-3. The re¬ 

combination rate, oB) is a function of temperature, which 
is calcula ted fr om the stellar spectrum and the clump 
density (' 112.2.31) . The free-free emission is evaluated by 
solving equation © assuming spherical geometry. The 
FWHM of the velocity distribution of the ionized clump 
gas, which will be used as part of the diagnostic process 
to classify H II regions (Q, is estimated simply from 

thermal broadening, ufwhm.cI = VS In 2cr„. This esti¬ 
mation gives the minimum value of the line width, since 
the clump gas could also be disturbed by outflows and 
additional turbulent motions that broaden the line pro¬ 
file. 

We note some caveats of this Stromgren sphere ap¬ 
proximation for ambient clump ionization. First, the 
geometry of the ionized region is likely to be affected 
by the outflow cavity, leading to deviations from a sim¬ 
ple sphere. Second, the ambient clump density is ex¬ 
pected to be inhomogeneous, which can lead to various 
effects, including enhanced effective densities at ioniza- 
tion fronts that can lead to smaller H H region sizes (e.g., 
iTan &: McKeel[200ll . Thus, the free-free flux evaluated 
by the Stromgren sphere analysis is only expected to be 
a very crude approximation. If the ambient clump region 
is optically thin, its total free-free flux from the ionized 
region is just proportional to the total escaping photon 
rate, /esc'S'*,acc, regardless of its shape or density. How¬ 
ever, if the ionized region is optically thick, the free-free 
flux is smaller than the optically thin limit by the factor 
of its mean optical depth. Therefore, the error in the 
free-free flux would be small if the ionized clump region 
is optically thin. 

4. DISCUSSION 

4.1. Outflow-confined HII regions and their break out 


Here we summarize the evolution of ionized outflows. 
Note, the increase of ionizing flux in the Kelvin-Helmholz 
contraction phase is very rapid (Fig. [2]) and our numer¬ 
ical results have only been carried out with relatively 
coarse sampling across the evolutionary sequence. How¬ 
ever, we can still present a description of the following 
evolutionary sequence. At the first stage, when the ion¬ 
izing flux is small, the MHD-driven outflow is the first 
barrier to the propagation of ionizing photons, except 
in the very narrow range of directions within the inner 
disk wind cavity. The H H region is confined in a very 
thin layer of the inner outflow wall, e.g., model Ahl6 in 
Figure [51 As the ionizing flux increases, the H II region 
extends to greater distance from the protostar. At a crit¬ 
ical ionizing flux, the H II region breaks out at a middle 
latitude (see model A12 in Fig|4]). At the lower latitude, 
the density is higher than that in the middle latitude 
since it is closer to the launching zone of the disk wind. 
On the other hand, the density gradient at the higher 
latitude is flatter than that at the middle latitude due 
to the collimated structure of the wind. Therefore, the 
H H region breakout starts from 6 ~ 6(w,esc/2. The crit¬ 
ical ionizing flux for the radial breakout depends on the 
detailed structure of outflow, i.e., mass loss rate, opening 
angle and velocity structure. From a qualitative stand¬ 
point, the mass loss rate is roughly proportional to the 
mass accretion rate, and thus the critical ionizing flux for 
the breakout is higher for higher-Sd cases. The breakout 
does not occur at 8Mq even in the low-Sci case of model 
A108. The typical flux for breakout is s“^ in 

our numerical calculations, which is as high as the flux 
from the ZAMS at about 10-15 Mq. All models in this 
work reach 20 Mq at the final stage, however, it should 
be noted that there are almost as many stars between 8 
and 10 Mq as above 20 Mq. In the formation of such 
stars with < 10 Mq, the H H region is confined during 
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Figure 12. The SEDs at inclination angles of 60° (left) and 30° (right) for model group A, Al, Ah, B and C (from top to bottom). The 
dotted lines show the free-free emission components. 


the main accretion phase and breakout occurs only to¬ 
wards the end of accretion when the outflow rate also 
decreases. 

After first breakout at 0w,esc/2, the photoionized re¬ 
gion extends in polar angle. To ionize almost the en¬ 
tire outflow, the ionizing flux needs to increase by about 
one order of magnitude from the first breakout ionizing 
flux. This takes only 10^-10'^ yr since the ionizing flux 
increases rapidly in the KH contraction phase. Once the 
entire outflow is ionized, the expansion in polar angle 
of the H II region slows down matching the increase of 
the outflow cavity opening angle, which has a timescale 
similar to the accretion time, i.e., m*/m* ^ 10^-10® yr. 

Our model of photoionized disk winds may 
also help to explain the short timescale variabil- 
ity that is seen in some UC/HC H II regions 
(iFranc o-Hern andez &: Rodriguez! 2004 IRodrfguez et al.l 
[20071 iGalvan-Madrid et al.lT200S ) and which has also 


been model ed in the r a diativ e hydrodynamical sim¬ 
ulations of iPeters et al.l (|2010D . Fluctuations in the 
accretion rate would also lead to changes in the outflow 
massloss rate, which would change the density and thus 
the extent of the photoionized region. Changes could 
happen on outflow crossing timescales, i.e., as short as 
^ 0.01 pc / 1000 km s~^ ~ 10 yr. 

Our model makes prediction for the spectral and the 
size indices, i.e., p and q, of the free-free emission. The 
free-free emission spectral index at 10 GHz is about 
p = 0.4-0.7 (see §3.2.1), which indicates that the out¬ 
flow is partially optically thick. The apparent size of 
the observed free-free emission (see §3.2.2) is as small as 
^ 10“^ pc at 10 GHz and smaller at higher frequencies, 
scaling with size as r oc u'J, with the index of q ~ —0.7, 
even though the ionized region is larger than the core 
scale of 0.1 pc after breakout. These frequency depen¬ 
dences are related to the density structure of the ionized 
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Figure 13. Resolved (top) and convolved (bottom) images of 
model A16 (Me = 60 Mq,Ec 1 = 1 gcm~^, and m* = 16 Mq) at 
an inclination angle of 60°. The panels show the images at the fre¬ 
quencies of 1.5, 10, and 870 GHz from left to right. In the top pan¬ 
els, the box sizes are 20, OOOAU x 20, OOOAU, 4, OOOAU x 4, OOOAU, 
and lOOAU x lOOAU, respectively. In the bottom panels, the beam 
sizes (FWHM) at each frequency are indicated at the bottom-left 
corners (1.3, 0.2, and 0.01"). A distance of 1 kpc is assumed. In the 
convolved images of the bottom row, from left to right the bright¬ 
ness temperature of is equivalent to 3.1/iJy(Tji,/K) beam~^, 
3.3pJy(Tb/K) beam“^ and 62pJy(Tb/K) beam“^ at these given 
frequencies. 

outflow. A simplest analytical structure is a spherical 
wind with constant velocity, i.e., riH = di„,/(47rr^z;^) 
where m,„ and v,„ are the w ind mass loss rate and ve¬ 
locity (|Panagia &: FelliiriQTSD . Since the opacity for the 
free-free emission at radio frequencies is approximatly 
oc if the gas is isothermal, the optical depth 

rapidly decreases as r~^. Therefore, the outer region be¬ 
comes optically thin and has less contribution to the to¬ 
tal free-free flux even though the ionized gas extends to 
larger scales. In this spherical model, the radius of the 
optically thick region scales as oc , and thus 

the total free-free flux, which is proportional to the flux 
from the optically thick re gion, is oc oc 

Additionally, as shown bv iRevnoldsj (|1986ll . even with a 
bipolar jet geometry, the set of indices will be similar to 
the simple spherical wind if the jet is conical (constant 
aperture angle), of constant velocity and ionization frac¬ 
tion, and is isothermal. Thus the departure from the 
indices of (p, q) = (0.6, —0.7) is related to acceleration, 
collimati on, thermal structure, or variation o f ionization 
fraction (|Revnoldslll986t [Guzman et al.ll20l4l . Our wind 
model is more complex compared to the spherical wind 
or the conical jet, however the fundamental physics is the 
same and thus the values of its spectral and size indices 
are similar. 

4.2. Comparison with observed radio sources around 
massive protostars 

Here, we compare general properties of our model with 
observed radio sources around massive protostars and 
show that the features of our theoretical models are sim¬ 


ilar to those of observed radio winds and jets. 

Figure ITSl shows the properties of size, line width, ra¬ 
dio and bolometric luminosities, and outflow momentum 
rate of our models and of observed radio sources. For 
our models, size and free-free luminosity are measured 
at 8 GHz and line width is evaluated from H40a assum¬ 
ing an inclination angle of 60°. The size is evaluated by 
the diameter at which half of the total radio flux of the 
system is emitted. As we have seen, the inclination angle 
does not affect these properties significantly at radio fre¬ 
quencies. Note, the outflow momentum rate, p^, of our 
model shown in Figure[T8]is obtained from the protostel- 
lar evolution calculation of ZTH14 rather than from an 
estimation via a synthetic observation. 

In p anel (a) of Figure [181 show how iHoare et al.l 
(l2007ri distinguished objects into three categories based 
on size and line width: UG H II regions are ~ 0.1 pc 
and < 40 kms“^, HG H II regions are < 0.05 pc and 
> 40 kms“^, and massive young stellar object (MYSO) 
wind sources are < 0.005 pc and > 70kms“^. It should 
be noted that this classification of the radio sources is 
somewhat arbitrary. For example, iKurtzl (I2002L 120051) 
defined UG H II regions as having diameters < 0.1 pc 
and HC H II regions as < 0.01 pc w ith emission mea¬ 
sures > 10® pc cm“®. iTan et ah! (|2014D adopts a simpler 
definition only based on the abov e siz es. When we re¬ 
fer to data from iHoare et al.l (|2007t) and IHoare fc Francol 
ll200l . we follow Hoare et al.’s classification. The line 
width of UG and HC H II regions are derived from radio 
HRLs, while those of wind sources are from IR HRLs. 
The wind sources (BN, W33A, S140 IRS 1, NGC 2024 
IRS 2, GL 989 and GL 490 are plotted here) are the 
smallest and fastest objects. Comparing in those plots, 
our models reproduce the properties of wind sources 
quantitatively, with sizes of ^ 10“® pc and line widths 
of ^ 100kms“^. Two of our models have broader 
line widths than 250kms“^ (718kms“^ for A12 and 
429kms“^ for C16). These two models are at the be¬ 
ginning of H II region expansion, and the ionized region 
is confined only near the axis where the wind velocity is 
highest. Thus, their line width is exceptionally broad. 

The difference between UC/HC H II regions and 
MYSO wind sources is more clearly seen in the ratio of 
radio and bolometric luminosities in panel (b) of Figure 
[T8l UC/HC H II regions are more radio-loud than MYSO 
wind sources by two orders of magnitude or more. Even 
though the average radio to bolometric luminosity ratio 
of our models are about one order of magnitude higher 
than those of MYSO wind sources, there is still overlap 
gi ven their dispe r sions . 

lAnglada et al.l (I2015D found two empirical correlations 
related to the radio luminosity of observed radio jets: one 
is a correlation with bolometric luminosity; the other is 
with outflow momentum rate: 


mjy kpe^ 


= 8 X I0“® 



F,d? 

mJy kpe^ 


= 190 


Pw 


Mq yr“^ km s 


-1 


0.9 


(17) 

(18) 


Interestingly, those two correlations are consistent over a 
wide bolometric luminosity range of ~ 1 Lq to ~ 10® Lq 
(dotted lines in panels (c) and (d) of Fig. [181). MYSO 
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Figure 14. Resolved (top) and convolved (bottom) images at 10 GHz of all models at inclination angles of line of sight of 60° (left) and 
30° (right) to the outflow axis. The box size of top panels are 4, 000 AU x 4, 000 AU. The images in the bottom panels are convolved with 
a beam size of 0.2" and a distance of 1 kpc is assumed. Here, the brightness temperature of Tf, is equivalent to 3.3/xJy(Tjij/K) beam”^. 


wind sources from iHoare fc Francol (|2007li (black cir¬ 
cles) also fall on the first correlation. The solid line 
in panel (b) is the expected radio luminosity in the 
optically thin limi t using the ioniz ing photon rate of 
ZAMS models by iThomosonl (jl984f l. which represents 
the maximum radio luminosity from photoionized gas. 
At Lboi 10^ Lq, however , there are radio jets that ex¬ 
ceed this maximum value (jAnglada et al.l[l992f) . It has 
been suggested that the single power law correlations 
of equations dm) and dD indicates a common nature 
of (ioniz ed) outflows from low-mass and high-m ass pro- 
tostars (jR.odriguez et al.lll9891j : lAngladal 1199,^ . Since 
the photoionization properties differ greatly from low- 


mass to high-mass stars, this might imply that a different 
mechanism, i.e., shock ionization, is responsible for the 
observed radio emission along these correlations. How¬ 
ever, we find that the radio luminosities of our model 
agree reasonably well with both the correlations of eqs. 
(Hzl and (|18|1 , even though they are photoionized struc¬ 
tures. We predict that there should be a transition to 
photoionization-dominated winds for massive protostars, 
but that in this regime the radio/bolometric luminosities 
and the outflow momentum rate may not differ too much 
from this empirical relation. The systematic differences 
of our models are towards somewhat higher values of ra¬ 
dio luminosities by factors of 1.5-13. This may be a 
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Figure 15. H40 q: line profiles of model A16 viewed with in¬ 

clination angles of 60° (red solid line) and 30° (blue solid line). 
The FWHM are ufwhm = 104, 115kms“^ at 60, 30° respec¬ 
tively. The gray dotted line shows a Gaussian distribution with 
upwHM = 110 kms“^ and a peak flux of 28 mjy. A distance of 
1 kpc is assumed. 


-125kni/s Okm/s 125 km/s 



logTb[K] 

Figure 16. H40a maps of model A16 at—125, 0, and 125km s“^ 
(from left to right). Top two rows show the resolved images viewed 
at inclination angles of 60° (first row) and 30° (second row). The 
box size is 1500AU x 1500AU. Bottom two rows show the convolved 
images viewed at inclination angles of 60° (third row) and 30° 
(fourth row). A beam size of 0.1'' and a distance of 1 kpc are 
assumed. 




Figure 17. Top: Escape fraction of ionizing photons to each 
polar angle for models A12, A16, and A24 (red, blue, and black 
solid lines, respectively). The dotted lines represent the opening 
angle of each model. Bottom: The total escape fraction of ionizing 
photon vs. the opening angle of the disk wind cavity. The circles 
show the models and the dotted lines represent their evolution 
(purple: model group A; green: Al; light blue: Ah; orange: B; and 
yellow: C). The size of each circle represents the protostellar mass. 
The solid line shows the covering fraction of the whole sky solid 
angle by the disk wind cavity, i.e., 1 — cos0w, esc. 

signature of the initial stage of photoionization leading 
to enhanced radio fluxes. Note, our models tend to be 
more massive than the observed sources, so more data in 
this range is needed. 

If the photoionized region is optically thin, the free-free 
emission is simply proportional to the ionizing photon 
flux rate ([Hummer fc SeatCTll963l) . However, our model 
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Figure 18. Properties of theoretical outflow-confined H I I region models c ompa red with observed radio sources around protostars. 
The observational da ta are adopted from (Hoare et aTI l|2007l) . (Hoare & Franco! (|2007l ) (filled circles in panels (a), (b), and (c)), and from 
lAnglada et al.l II2015I) (filled triangles in panels (b) and (c)). In all panels, the theoretical models are shown by square symbols with sizes 
indicating the stellar mass m* (purple: model group A, green: Al, light blue: Ah, orange: B, and yellow: C). The dash lines represent 
evolution within the model groups, (a): Size versus line width. Observed objects are shown by circles (blue: UC H II regions, red: HC 
H II regions, and black: MYSO wind sources). The sizes for the theoretical models (which vary with frequency) and the observed sources 
have both been evaluated at 8GHz. (b): The ratio of radio and bolometric luminosities versus line width. Symbols are same as the panel (a). 
The radio luminosity for the theoretical models, which varies with frequency, and the observed sources have been evaluated at 8 GHz. (c): 
Radio luminosity versus bolometric luminosity. Triangles show observed radio jets (light green: high-luminosity, and blue: low-luminosity 
sources). The radio luminosities of the theoretical models are evaluated at 8GHz, and the observational data are at 5 or 8GHz. The dotted 
line shows the fitting by lAnglada et all II2015I) : Thd^/mJy kpc^ = 8 x 10~^ (Lbol/^©)^'^- This relation is based on observed radio jets over 
a wide bolometric luminosity range of 1-10^ Lq. Th e solid line shows the estimated radio flux from an optically thin H II region based on 
the ionizing luminosities of ZAMS stellar models by IThompsonI (119841) . The purple line shows the estimated radio flux from an optically 
thin H II region based on the ionizing luminosities of the protostar of model group A. (d): Radio luminosity versus outflow momentum 

rate. The dotted line shows the fitting by lAnglada et all J2015I) : Ft/d^/mJy kpc^ = 190 (pw/M q yr“^ km The square and triangle 

symbols and the purple line are same as in panel (c). The momentum rate of our model pw is obtained from the evolution calculation by 
ZTH14 rather than from estimation by a synthetic observation. Properties of the theoretical models of outflow-confined H II regions are 
broadly consistent with those of the high luminosity end of the observed radio wind/jet sources. 
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indicates much lower free-free flux than the optically thin 
limit (panel (c) in Fig. ITSl) . First, this is because more 
than half of the emitted ionizing photons are absorbed 
at the innermost region of r <C lOOOAU. This region is 
optically thick for the radio free-free emission, and thus 
the absorption there does not contribute to the total free- 
free flux. In the case of model A16, 78% of photons 
are absorbed inside r < lOOAU. Second, most of the 
rest of the photons, /esc, escape from the system. Those 
escaping photons do not induce free-free flux from the 
core scale, but contribute to the flux from the ambient 
clump gas, which is discussed below 1 114.31) . The final 
reason is that the photoevaporation flow has so far not 
been included in our model. If the ionizing photon rate 
is high enough to ionize the disk and envelope surfaces, 
e.g., in the case of model A24, the outflow rate of ionized 
gas increases and the free-free flux should be enhanced. 
While the first and second reasons are physical, the last 
one is due to the incompleteness of our model. We will 
address photoevaporation and its feedback effects in a 
subsequent paper. 

Our ionized outflow models typically have spectral in¬ 
dices with p = 0.4-0.7 (Fig. [T21), and size index with 
frequency of g ~ —0.7 from free-free emission at 10 GHz 
(Fig. [T31). As discussed above, those indices are related 
to the density gradient of the outflow. The indices ob¬ 
tained by our calculations are consistent with some ob¬ 
served radio jets, such as MWC 349 (p = 0.67 ± 0.03 and 
q = -n.74-kn.03. lTafova et al.][2nnl and W75N(B') (p = 
0.61T 0.03 and q = —0.7±0.03. lCarrasco-Gonzaiez et aP 
I20I5I1 . G345.4938-1-01.4677, however, has the slightly 

steeper indices of p = 0.92 ± 0.01 and q = —I.IIqI 
than our model, which may be interpreted as being 
due to a stronger deg ree of collimation and acceleration 
(jGuzman et al.l[^14 Guzman et ah, in prep.). 


We should notice that, even though our theoretical 
model has similar properties to those of the observed 
radio jets, there is still a lack of some physical pro¬ 
cesses in our model. The velocity of some jets is esti¬ 
mated to be as high as 500 kms“^ from the proper mo- 
tion of associated l obes along the axes, e.g., HH 80-81 
(|Marti et al.lll998L 119951: fMasoue et al.1 120151) and Cep 
A HW 2 (jGuriel et al.ll2006l l. Similar to these sources, 
our wind models have high radio brightnesses and high 
velocities of > 300 kms“^ near the axes (Figs. [1] [71 
etc.). However, unlike our model, some of lobes are 
found to have the radio e mission with negative sp ectral 
indices , such as Serpens (|Bodriguez et al.l 1989all. HH 
80-81 (|Marti et al.l 1199311. Cenh A H W 2 ( Garav et al.l 
[T9M1 . W3(OHl (IWilner eir^HO^ . and IRAS 16547- 
4247 ([Rodriguez et al.l 1200511 . These negative spectral 
indices can be interpreted as due to nonthermal syn¬ 
chrotron e mission induced by s trong shocks with the am¬ 
bient gas ([Araudo et al.l l20?i^ . and this has been con¬ 
firmed by the detection of linearly polari zed emission in 
HH 80-81 (iCarrasco-Gonzalez et al.ll2010ll . Additionally, 
our simple axisymmetric model has not accounted for the 
possibility of outflow axis precession, which is apparent 
in IRAS16547-4247 ([Rodriguez et al'l [200811 . Although 
the current number of sources with these features is rel¬ 
atively small, inclusion of such additional physics in our 
theoretical models, i.e., shock ionization and outflow axis 
precession, are likely to be needed for a complete under¬ 


standing of jet properties in massive star formation. Such 
analysis is also deferred to a future paper. 

As we have seen, the properties of our models are 
con sistent with the so -called “MYSO wind s ources ” 
in iHoare et al.l (l2007| l and iHoare fc Francol ([200711 . 
and “radio je t s” wi th high bolometric luminosity in 
lAnglada et al.l (1201511 . In general, the main difference 
of “winds” and “jets” is the shape, i.e., winds have rela¬ 
tively wide opening angles and jets are more collimated 
bipolar outflows. This difference could be related to evo¬ 
lutionary stage (see also Fig. dl. To clarify our view of 
the evolutionary sequence of massive star formation and 
its outflow launching mechanism, further study is needed 
utilizing multi-wavelength observations to be compared 
with theoretical models. 

4.3. Ionization of ambient clump gas by escaping 
photons 

Finally, we comment on the ionization of the ambient 
clump. Tens of percent of the total ionizing photon out¬ 
put escape along the outflow cavity (Fig. [IT]). These 
photons are likely to ionize the ambient clump which 
would lead to enhanced radio luminosity. We estimate 
fea tures of these extended ionized regions, as described 
in 93.41 Figure |T9| shows the comparison of our ionized 
clump models and observed sources (the observational 
data are the same as in Fig. m- The Stromgren-based 
models of ionized ambient clumps typically have the size 
of ^ 10“^ pc and line widths of ^ 20kms“^, and large 
radio luminosities of > 10^ mJy kpc^, which are larger in 
size, narrower in line width, and more radio loud than 
the ionized outflow models. These features are roughly 
consistent with observed UG/HG H II regions, except the 
line widths are broader in the observed systems. How¬ 
ever, this difference is not surprising as the theoretical 
models include only thermal broadening and so give min¬ 
imum widths. 

We conclude that some UC/HC H II regions may be 
formed by the escaping photons from massive protostars 
with radio jets or winds. In this study, the ionized clump 
features are evaluated based on the simple Stromgren 
analysis ( 93.41) . For more detailed analysis, planned for 
future work, one needs to model more accurately the 
structure and kinematics of ambient clump gas, including 
its perturbation by the outflow, to better understand the 
connection of radio winds/jets and UG/HG H II regions. 

5. SUMMARY 

We have studied the structure, evolution and observa¬ 
tional properties of photoionized disk wind outflows that 
are expected to be associated with massive protostars 
forming by Gore Accretion. Specifically, we have followed 
the quantitat ive properties o f the T urbulent Core Accre¬ 
tion Model of lMcKee fc TanI ([20031 ) and follow-up proto- 
stellar evolution, disk wind and IR con tinu um radiative 
transfer mod els of I Zhang fc TanI (I2011D and IZhang et al.l 

dMllloil). 

We calculated the transfer of ionizing photons and ob¬ 
tained the photoionized structure. The disk wind starts 
to be ionized when the stellar mass reaches about 10 - 
20 Mq , depending on the properties of the initial cloud 
core, due to the stellar temperature increase resulting 
from Kelvin-Helmholz contraction. At this moment, the 
photoionized region is confined to zenith angles of about 
























































Outflow-Confined HII Regions 


21 


[a] [b] 




FWHM [km/s] FWHM [km/s] 


[C] 



log Lboi [Lq] 


Figure 19. Same as panels (a), (b) and (c) of Figure ITsl except 
the disk wind cavity. 

10 - 30°, while the high density regions near the outflow 
axis and equatorial plane remain neutral. On the other 
hand, in the radial direction, the photoionized region is 
unconfined, since the density gradient of the disk wind is 
steep. The disk wind is almost entirely ionized in about 
10^ - 10“^ yr, and the photoionized region extends as the 
outflow opening angle increases. The outer disk surface 
is not ionized even though the diffuse radiation has been 
considered. This is because of the shielding effects of the 
inner disk and the disk wind. 

Using the obtained photoionized structures, we calcu- 


theoretical models now account for the escape of ionizing photons from 

lated radiative transfer of free-free continuum and Hydro¬ 
gen recombination line emission. The free-free emission 
exceeds dust emission at the radio frequencies lower than 
100 GHz once the ionized region is formed. While shorter 
wavelength infrared flux from dust emission can depend 
on inclination of viewing angle, the radio flux from free- 
free emission is almost independent of inclination and its 
typical flux density is ^ (20 - 200) x (u/lOGHz)^’mJy 
at a distance of 1 kpc with a spectral index p ~ 0.4 - 
0.7 The apparent size depends on frequency as approxi¬ 
mately 500(u/10GHz)“°-^ AU. The profile of H40a has 
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a FWHM of ^ 100 km with broad wings. These 
properties of our model are similar to those of observed 
radio winds and jets from massive protostars. Detailed 
comparison with individual sources is deferred to future 
papers in this series. 

ACKNOWLEDGMENTS 

The authors thank Christopher McKee, Andres 
Guzman, Taishi Nakamoto, Hiroyuki Kurokawa, Kohei 
Inayoshi, and Shuo Kong for fruitful discussions and com¬ 
ments. The authors also thank the anonymous referee 
for comments, which were useful to improve the original 
manuscript. JCT acknowledges support from NSF grant 
AST 1411527. 

REFERENCES 


Anglada, G., Rodriguez, L. F., Canto, J., Estalella, R., Sc 
Torrelles, J. M. 1992, ApJ, 395, 494 
Anglada, G. 1995, Revista Mexicana de Astronomia y Astrofisica 
Conference Series, 1, 67 

Anglada, G., Rodriguez, L. F., S>c Carrasco-Gonzalez, C. 2015, 
Advancing Astrophysics with the Square Kilometre Array 
(AASKA14), 121 

Araudo, A. T., Romero, G. E., Bosch-Ramon, V., Sz Paredes, 

J. M. 2007, A&A, 476, 1289 

Baez-Rubio, A., Martin-Pintado, J., Thum, C., S>c Planesas, P. 
2013, A&A, 553, AA45 

Beuther, H., Schilke, P., Sridharan, T. K., et al. 2002, A&A, 383, 
892 

Blandford, R. D., & Payne, D. G. 1982, MNRAS, 199, 883 
Bunn, J. C., Hoare, M. G., Sc Drew, J. E. 1995, MNRAS, 272, 346 
Butler, M. J., Sc Tan, J. C. 2012, ApJ, 754, 5 
Butler, M. J., Tan, J. C., Sc Kainulainen, J. 2014, ApJ, 782, LL30 
Carrasco-Gonzalez, C., Rodriguez, L. F., Anglada, G., et al. 2010, 
Science, 330, 1209 

Carrasco-Gonzalez, C., Torrelles, J. M., Canto, J., et al. 2015, 
Science, 348, 114 

Castelli, F., Sc Kurucz, R. L. 2004, lAU Symp. No 210, Modelling 
of Stellar Atmospheres, eds. N. Piskunov et al. 2003, poster 
A20, arXiv:astro-ph/0405087 

Codella, C., Cabrit, S., Gueth, F., et al. 2014, A&A, 568, L5 
Curiel, S., Ho, P. T. P., Patel, N. A., et al. 2006, ApJ, 638, 878 
Dale, J. E., Bonnell, 1. A., Clarke, C. J., Sc Bate, M. R. 2005, 
MNRAS, 358, 291 

Davies, B., Hoare, M. G., Lumsden, S. L., et al. 2011, MNRAS, 
416, 972 

Draine, B. T. 2011, Physics of the Interstellar and Intergalactic 
Medium by Bruce T. Draine. Princeton University Press, 

2011. ISBN: 978-0-691-12214-4 

Dupree, A. K., Sc Goldberg, L. 1970, ARA&A, 8, 231 
Ferland, G. J., Porter, R. L., van Hoof, P. A. M., et al. 2013, 
RMxAA, 49, 137 

Figer, D. F. 2005, Nature, 434, 192 

Franco-Hernandez, R., Sc Rodriguez, L. F. 2004, ApJ, 604, L105 
Galvan-Madrid, R., Rodriguez, L. F., Ho, P. T. P., Sc Keto, E. 
2008, ApJ, 674, L33 

Garay, G., Mardones, D., Bronfman, L., et al. 2007, A&A, 463, 
217 

Garay, G., Ramirez, S., Rodriguez, L. F., Curiel, S., Sc Torrelles, 
J. M. 1996, ApJ, 459, 193 

Gibb, A. G., Hoare, M. G., Little, L. T., Sc Wright, M. C. H. 

2003, MNRAS, 339, 101 

Goodman, A. A., Benson, P. J., Fuller, G. A., Sc Myers, P. C. 
1993, ApJ, 406, 528 

Guzman, A. E., Garay, G., Brooks, K. J., Sc Voronkov, M. A. 

2012, ApJ, 753, 51 

Guzman, A. E., Garay, G., Rodriguez, L. F., et al. 2014, ApJ, 

796, 117 

Hoare, M. G., Kurtz, S. E., Lizano, S., Keto, E., Sc Hofner, P. 
2007, Protostars and Planets V, 181 


Hoare, M. G., Sc Franco, J. 2007, A Volume Honouring John 
Dyson, Edited by T.W. Hartquist, J. M. Pittard, and S. A. E. 
G. Falle. Series: Astrophysics and Space Science Proceedings. 
Springer Dordrecht, 2007, p.61, arXiv:0711.4912 
Hollenbach, D., Johnstone, D., Lizano, S., Sc Shu, F. 1994, ApJ, 
428, 654 

Hosokawa, T., Sc Omukai, K. 2009, ApJ, 703, 1810 
Hosokawa, T., Yorke, H. W., & Omul^i, K. 2010, ApJ, 721, 478 
Hosokawa, T., Omukai, K., Yoshida, N., & Yorke, H. W. 2011, 
Science, 334, 1250 

Hummer, D. G., Sc Seaton, M. J. 1963, MNRAS, 125, 437 
Keto, E. 2007, ApJ, 666, 976 

Klaassen, P. D., Juhasz, A., Mathews, G. S., et al. 2013, A&A, 
555, A73 

Konigl, A., Sc Pudritz, R. E. 2000, Protostars and Planets IV, 759 
Kratter, K. M., Matzner, C. D., Sc Krumholz, M. R. 2008, ApJ, 
681, 375 

Krumholz, M. R., Klein, R. L, McKee, C. F., Offner, S. S. R., Sc 
Cunningham, A. J. 2009, Science, 323, 754 
Krumholz, M. R., McKee, C. F., Sc Klein, R. 1. 2005, ApJ, 618, 
L33 

Krumholz, M. R., Stone, J. M., Sc Gardiner, T. A. 2007, ApJ, 

671, 518 

Kuiper, R., Klahr, H., Beuther, H., Sc Henning, T. 2010, ApJ, 

722, 1556 

Kuiper, R., Klahr, H., Beuther, H., Sc Henning, T. 2011, ApJ, 

732, 20 

Kurtz, S. 2002, Hot Star Workshop HI: The Earliest Phases of 
Massive Star Birth, 267, 81 

Kurtz, S. 2005, Massive Star Birth: A Crossroads of Astrophysics, 
227, 111 

Lanz, T., Sc Hubeny, L 2007, ApJS, 169, 83 
Li, J., Wang, J., Gu, Q., Zhang, Z.-y., Sc Zheng, X. 2012, ApJ, 
745, 47 

Lopez-Sepulcre, A., Codella, C., Cesaroni, R., Marcelino, N., Sc 
Walmsley, C. M. 2009, A&A, 499, 811 
Marti, J., Rodriguez, L. F., Sc Reipurth, B. 1993, ApJ, 416, 208 
Marti, J., Rodriguez, Sc Reipurth, B. 1995, ApJ, 449, 184 
Marti, J., Rodriguez, L. F., Sc Reipurth, B. 1998, ApJ, 502, 337 
Martins, F., Sc Plez, B. 2006, AScA, 457, 637 
Masque, J. M., Rodriguez, L. F., Araudo, A., et al. 2015, ApJ, in 
press (arXiv: 1510.01769) 

Matzner, C. D., & McKee, C. F. 2000, ApJ, 545, 364 
McKee, C. F., Sc Tan, J. C. 2003, ApJ, 585, 850 (MT03) 

McKee, C. F., Sc Tan, J. C. 2008, ApJ, 681, 771 
McLaughlin, D. E., Sc Pudritz, R. E. 1997, ApJ, 476, 750 
Meynet, G., Georgy, C., Revaz, Y., et al. 2010, Revista Mexicana 
de Astronomia y Astrofisica Conference Series, 38, 113 
Nakano, T., Hasegawa, T., Sc Norman, C. 1995, ApJ, 450, 183 
Nisini, B., Smith, H. A., Fischer, J., Sc Geballe, T. R. 1994, A&A, 
290, 463 

Palau, A., Puente, A., Girart, J. M., et al. 2013, ApJ, 762, 120 
Panagia, N., Sc Felli, M. 1975, A&A, 39, 1 

Peters, T., Klessen, R. S., Mac Low, M.-M., Sc Banerjee, R. 2010, 
ApJ, 725, 134 

Peters, T., Banerjee, R., Klessen, R. S., & Mac Low, M.-M. 2011, 
ApJ, 729, 72 

Reynolds, S. P. 1986, ApJ, 304, 713 
Richling, S., Sc Yorke, H. W. 1997, A&A, 327, 317 
Rodriguez, L. F., Curiel, S., Moran, J. M., et al. 1989a, ApJ, 346, 
L85 

Rodriguez, L. F., Garay, G., Brooks, K. J., Sc Mardones, D. 2005, 
ApJ, 626, 953 

Rodriguez, L. F., Gomez, Y., Sc Tafoya, D. 2007, ApJ, 663, 1083 
Rodriguez, L. F., Moran, J. M., Franco-Hernandez, R., et al. 

2008, AJ, 135, 2370 

Rodriguez, L. F., Myers, P. C., Cruz-Gonzalez, L, Sc Terebey, S. 
1989b, ApJ, 347, 461 

Salem, M., Sc Brocklehurst, M. 1979, ApJS, 39, 633 
Shakura, N. L, Sc Sunyaev, R. A. 1973, A&;A, 24, 337 
Shu, F. H. 1977, ApJ, 214, 488 

Stone, J. M., Mihalas, D., Sc Norman, M. L. 1992, ApJS, 80, 819 
Strelnitski, V. S., Ponomarev, V. O., Sc Smith, H. A. 1996, ApJ, 
470, 1118 

Tafoya, D., Gomez, Y., Sc Rodriguez, L. F. 2004, ApJ, 610, 827 


Outflow-Confined HII Regions 


23 


Tan, J. C., & McKee, C. F. 2003, lAU Symposium 221, Star 
Formation at High Angular Resolution, Eds. M. Burton, R. 
Jayawardhana, Sz T. Bourke, arXiv:astro-ph/0309139 
Tan, J. C., & McKee, C. F. 2004, ApJ, 603, 383 
Tan, J. C., & McKee, C. F. 2004, The Formation and Evolution 
of Massive Young Star Clusters, 322, 263 
Tan, J. C., Beltran, M. T., Caselli, P., et al. 2014, Protostars and 
Planets VI, 149 

Tanaka, K. E. L, Sz Nakamoto, T. 2011, ApJ, 739, L50 
Tanaka, K. E. L, Nakamoto, T., Omukai, K. 2013, ApJ, 773, 
155 

Thompson, R. I. 1984, ApJ, 283, 165 

Thum, C., Matthews, H. E., Martin-Pintado, J., et al. 1994, 
A&A, 283, 582 


Ulrich, R. K. 1976, ApJ, 210, 377 
Walmsley, C. M. 1990, A&AS, 82, 201 
Weingartner, J. C., Sz Draine, B. T. 2001, ApJ, 548, 296 
Wilner, D. J., Reid, M. J., Sz Menten, K. M. 1999, ApJ, 513, 775 
Wood, D. O. S., Sz Churchwell, E. 1989, ApJS, 69, 831 
Yorke, H. W., Sz Bodenheimer, P. 1999, ApJ, 525, 330 
Yorke, H. W., Sz Welz, A. 1996, A&;A, 315, 555 
Zhang, Y., Sc Tan, J. C. 2011, ApJ, 733, 55 
Zhang, Y., Tan, J. C., Sc McKee, C. F. 2013, ApJ, 766, 86 
Zhang, Y., Tan, J. C., Sz Hosokawa, T. 2014, ApJ, 788, 166 
(ZTH14) 


